function IdSet=OneSide4a(PYX_cond, u0grid, u1grid, u2grid,...
                        equalorder, ...
                        Ugridreduced, sgreduced,...
                        id_Ugridreduced,...
                        s_id_gridreduced, sU_id_gridreduced,...
                        n_U_sets, Tcoord, indices_pairs,...
                        numeqconstraints,coordrowseq, fillAeq)
%% Complete equality constraints for parts invariant across parameter values
beq=[PYX_cond(1)-1;... %observational equivalence
      PYX_cond(2)-1;...
      PYX_cond(3)-1;...
      zeros(9,1);... %zeros
      1]; %ones   %[#equalityconstraints]x1 

%% Loop over parameter values 
exists=zeros(sgreduced,1);

for g=1:sgreduced
% Set indices to keep track of equal elements     
idU1=id_Ugridreduced{1}(g,:);
idU2=id_Ugridreduced{2}(g,:);
idU1_t=id_Ugridreduced{1}(g,:).';
idU2_t=id_Ugridreduced{2}(g,:).';
% Relevant sizes of the vector of unknowns
s1=s_id_gridreduced(g,1);
s2=s_id_gridreduced(g,2);
sU=sU_id_gridreduced(g); 

%%%%%%%%%%%%%%%%%%%%%% Equality constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%         
coordcolumnseq=[changecoord(s1,s2,idU1(1),idU2(1),1) changecoord(s1,s2,idU1(4),idU2(1),1) changecoord(s1,s2,idU1(1),idU2(4),1) ...%observational equivalence
                changecoord(s1,s2,idU1(2),idU2(2),1) changecoord(s1,s2,idU1(4),idU2(2),1) changecoord(s1,s2,idU1(2),idU2(4),1) ...
                changecoord(s1,s2,idU1(3),idU2(3),1) changecoord(s1,s2,idU1(4),idU2(3),1) changecoord(s1,s2,idU1(3),idU2(4),1) ... 
                changecoord(s1,s2,idU1((1:5)),idU2(5),1) ... %zeros
                changecoord(s1,s2,idU1(5),idU2((1:4)),1) ...
                changecoord(s1,s2,idU1(4),idU2(4),1)]; %ones    
                
Aeq=sparse(coordrowseq, coordcolumnseq, fillAeq, numeqconstraints, sU);   %[#equalityconstraints]x[sU]

%%%%%%%%%%%%%%%%%%%%%%  Inequality constraints %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Repeat the way we constructed Tcoord_specific but now using the parameter values under consideration
vectors = cellfun(@(x) {x(g,:)}, Ugridreduced); %cell 1xn_U_sets
%For j=1,...,n_U_sets, it picks the g-th row of Ugridreduced{j}
T_temp = cell(1,n_U_sets);
[T_temp{:}] = ndgrid(vectors{:}); 
T_temp = cat(n_U_sets+1, T_temp{:});
T = reshape(T_temp,[],n_U_sets); %matrix of dimension as Tcoord reporting the elements of all the possible n_U_sets-tuples from the U-sets

D=[T(indices_pairs(:,1),:) T(indices_pairs(:,2),:) ...
   Tcoord(indices_pairs(:,1),:) Tcoord(indices_pairs(:,2),:)]; %n_U_sets-tuples n_U_sets-tuples coordinates coordinates


%Save in D1 the rows of D where the first duplet is <= then the second
%Only the columns 5:end (reporting column-coordinates) are relevant
D1=D(sum(D(:,1:2)<=D(:,3:4),2)==2, 5:end);
sd1=size(D1,1);
%Remember <= means: (<,<) or (<,=) or (=,<) or (=,=)

%Save in D2 the rows of D where the first duplet is >= (excluding =,=) then the second
%Only the columns 5:end (reporting column-coordinates) are relevant
D2=D((sum(D(:,1:2)>D(:,3:4),2)==2) |... %>,>
       (sum(D(:,1:2)>D(:,3:4),2)==1 & sum(D(:,1:2)==D(:,3:4),2)==1) , 5:end); %>,= or =,>
sd2=size(D2,1);

coordrowsineq=[kron((1:1:sd1), ones(1,4))  kron((sd1+1:1:sd1+sd2), ones(1,4))]; %1x[(sd1+sd2)*4]
%4 comes from the number of all possible vertices generated by each pair of duplets

coordcolumns1=[changecoord(s1,s2,idU1_t(D1(:,1)),idU2_t(D1(:,2)),1)  changecoord(s1,s2,idU1_t(D1(:,3)),idU2_t(D1(:,2)),1) ...
               changecoord(s1,s2,idU1_t(D1(:,1)),idU2_t(D1(:,4)),1)  changecoord(s1,s2,idU1_t(D1(:,3)),idU2_t(D1(:,4)),1)]; %sd1x4

coordcolumns2=[changecoord(s1,s2,idU1_t(D2(:,1)),idU2_t(D2(:,2)),1)  changecoord(s1,s2,idU1_t(D2(:,3)),idU2_t(D2(:,2)),1) ...
               changecoord(s1,s2,idU1_t(D2(:,1)),idU2_t(D2(:,4)),1)  changecoord(s1,s2,idU1_t(D2(:,3)),idU2_t(D2(:,4)),1)]; %sd2x4         
          
coordcolumnsineq=[reshape(coordcolumns1.', 1, 4*sd1) reshape(coordcolumns2.', 1, 4*sd2)];   %transform coordcolumns1,2 in a row vector by reshaping row-wise  
%1x[(sd1+sd2)*4]


fillAineq=[repmat([-1 1 1 -1], 1, sd1) repmat([-1 1 1 -1], 1, sd2)]; %remember: signs are flipped because inequalities are <= 
%1x[(sd1+sd2)*4]

bineq=zeros(sd1+sd2, 1);

Aineq=sparse(coordrowsineq, coordcolumnsineq, fillAineq, sd1+sd2, sU);   %[#inequalityconstraints]x[sU]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Run LP %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Gurobi
%model.A=[Aineq; Aeq];
%model.obj=zeros(sU,1);
%model.sense=[repmat('<', size(Aineq,1),1); repmat('=', size(Aeq,1),1)];
%model.rhs=full([bineq; beq]); 
%model.ub=ones(sU,1);
%model.lb=zeros(sU,1);
%params.outputflag = 0; 
%result=gurobi(model, params);
%if isfield(result,'x')
%    exists(g)=1;
%end   

% MOSEK
param_MOSEK.MSK_IPAR_LOG = 0; 
prob.blx=zeros(sU,1); %lower bound unknowns
prob.ulx=ones(sU,1); %upper bound unknowns
prob.c=zeros(sU,1); %objective function
prob.a=[Aeq; Aineq];
prob.blc=[beq; -Inf(size(Aineq,1),1)];
prob.buc=[beq; bineq];
[~,result]     = mosekopt('minimize echo(0)',prob, param_MOSEK);
if strcmp(result.sol.itr.prosta, 'PRIMAL_AND_DUAL_FEASIBLE') %if the primal is feasible
   exists(g)=1; 
end

end

%Enlarge results to the original grid of parameter vectors
exists=exists(equalorder);
IdSet = [u0grid(logical(exists),:)  u1grid(logical(exists),:)  u2grid(logical(exists),:)];

